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The crystal-melt interfaces of a binary hard-sphere fluid mixture in coexistence with a single- 
component hard-sphere crystal is investigated using molecular-dynamics simulation. In the system 
under study, the fluid phase consists of a two-component mixture of hard spheres of differing size, 
with a size ratio a — 0.414. At low pressures this fluid coexists with a pure fee crystal of the larger 
particles in which the small particles are immiscible. For two interfacial orientations, [100] and [111], 
the structure and dynamics within the interfacial region is studied and compared with previous 
simulations on single component hard-sphere interfaces. Among a variety of novel properties, it is 
observed that as the interface is traversed from fluid to crystal the diffusion constant of the larger 
particle vanishes before that of the small particle defining a region of the interface where the large 
particles are frozen in their crystal lattice, but the small particles exhibit significant mobility. This 
behavior was not seen in previous binary hard-sphere interface simulations with less asymmetric 
diameters. 



I. INTRODUCTION 

A fundamental understanding of the nucleation, 
growth kinetics and morphology of crystals grown from 
the melt requires a detailed -microscopic description of 
the crystal-melt interfacenuotl However, such interfaces 
are very difficult to probe experimentally and reliable ex- 
perimental data, especially for structure and transport 
properties, is rare. It is then not surprising that com- 
puter simulations have, in recent years, played a leading 
role in the determination of the microscopic structure, 
dynamics and thermodynamics of such systemsa. 

To date, the vast majority of simulation studies have 
focused on single component interfacial systems. Such 
studies rWHfte from simple model ^systems such as hard 
spheresBQolj or LennarcUJoap l 3c3 to WW "realistic" 
systems, r-juch as waterE3E3o, silicontSEj or simple 
metalal3'tl. In contrast, thexe-have been but few studies 
on multicomponent system d 19 H 2C l in spite of the fact that 
most materials of technological interest are mixtures (for 
example, doped semiconductors, alloys and intermetallic 
compounds). In such systems, the crystal and coexist- 
ing fluid have differing composition, in general, and the 
change in concentration as one traverses the interface 
from one bulk phase into the other becomes an object of 
study. 

Of particular interest to materials scientists is the 
degree of interfacial segregation - the preferential ad- 
sorption of one component (usually the "solute" ) at the 
interface. In addition, the phase diagrams for multi- 
component systems are significantly more varied and 
complex than single component systems due to the addi- 
tional dimension of concentration. For a binary system 
several types of solid-liquid equilibria are possible. If the 
two types of particles are similar, then one typically has 
coexistence between a binary fluid and a substitution- 
aly disordered solid of similar structure to that of the 



pure components. However, if the two types of parti- 
cles are substantially different in nature, then generally 
the binary fluid will either be immiscible in the pure 
coexisting solid, or will coexist with one or more or- 
dered crystal mixtures (e.g. intermetallic compounds). 
Previous simulation studies on binary crystal-melt inter- 
faces have exclusively focused on the former case, namely 
the equilibrium between thc_fluid and a disordered crys- 
tal. Davidchack and LairdtJ recently reported results 
for a binary hard sphere system in which a substitution- 
ally disordered face-centered-cubic (fee) crystal coexists 
with a binary fluid mixture. In a related study, Hoyt 
et al. examined the crystal- melt interface of a Cu/Ni 
mixturdla. In both studies the degree of solute segrega- 
tion was found to be negligible. 

In the two studies above, the disordered fee crystal 
was stabilized by the fact that the two components were 
quite similar in size -for example, in the hard-sphere sys- 
tem studied by Davidchack and Laird, the diameters of 
the two types of spheres making up the system differed 
only by 10%. In this work, we extend the previous stud- 
ies to hard-sphere mixtures with significant size asym- 
metry. For such systems, in which the diameters differ 
by more than about 85%, the disordered fee phase is 
no longer stable and only coexistence of the fluid with 
ordered crystal structures is possible. In this work we ex- 
amine the interface between a binary hard-sphere fluid 
mixture and a coexisting fee crystal in which the small 
particle is immiscible. 

Our system of choice is a binary hard-sphere mixture 
in which the ratio of the smaller particle diameter to 
that of the larger particle is 0.414. Hard spheres are 
an important reference system for the crystal-melt inter- 
faces of simple systems since the structure, dynamics and 
phase behavior of dense atomic systems are dominated 
by packing considerations with only minor influence from 
the attractive parts of the interactions. For example, 



2 



it has been recently demonstrated!^ that the interfacial 
free energy of close-packed metals can be quantitatively 
described using a purely hard sphere model. The specific 
diameter ratio of 0.414 was chosen because, to perform 
an interface simulation, accurate phase coexistence pa- 
rameters are required a priori, and the phase diagram for 
this binary system has been worked out via simulation 
in some detailed This phase diagram shows that at low 
pressures the fluid mixture coexists with a pure fee crys- 
tal of the larger particles, but that at higher pressures the 
crystal structure in equilibrium with the fluid is a 1:1 (or 
AB) "intermetallic" compound with an "NaCl" struc- 
ture (the small and large particles form interpenetrating 
fee lattices). (The existence of the "NaCl" structure at 
this diameter ratio had been predicted earlier, using cell 
theorynJ. The diameter ratio, a = 0.414 is necessary 
for an "NaCl" structure to attain its maximum packing 
fraction of 0.793.) Thus, this system allows us to study 
the interfaces between binary fluids and two types of or- 
dered crystal phases: single component and "NaCl" . In 
this work we present results for the former, but simu- 
lations on the fluid/ "NaCl" are under way and will be 
reported later. 



As mentioned in the introduction we are interested in 
the present study in the interface between an fee crystal 
consisting of pure large (type A) spheres and its coex- 
isting binary fluid at a diameter ratio, a = 0.414. The 
pressure-composition phase diagram for a binary hard 
sphere system with this diameter ratio has been previ- 
ously determined by Trizac and coworkersE3 and is shown 
in Fig. 
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II. DESCRIPTION OF THE SYSTEM 

We consider a two-component system consisting of 
hard spheres of differing diameters, given by a a and as- 
Without loss of generality, it is assumed that a a > ctb- 
The interaction between two particles of type i and j, 
€ {A, B}), respectively, is then given by 



<t>ij ( r ) 







r < Ui 
r > (7i 



(1) 



where r is the distance between the centers of the two 
interacting spheres, and cry is the distance of closest pos- 
sible approach. In addition, we define the spheres to be 
additive, that is, = (a, +<t,-)/2. The state of the sys- 
tem is then completely described by specifying the total 
density, p = pa + Pb = N/V, the mole fraction, xa, of 
the larger species, and the diameter ratio a = as /a a- 
Note that so defined one has a G (0, 1). In a single com- 
ponent system composed of hard spheres of diameter a 
the packing fraction, 77 (the fraction of the total volume 
occupied by the spheres) is given by, 
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where p is bulk density. For the binary hard sphere sys- 
tem described above the packing fraction is 
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FIG. 1: Pressure-concentration phase diagram of a binary 
hard-sphere system with a = 0.414. [Reprinted from Ref. 
22, by permission of the publisher, Taylor and Francis Ltd. 
(www.tandf.co.uk/journals).] Note that to make the phase 
coexistence lines easier to distinguish, the pressure is plotted 
against and not xa as in the usual case. 

For this study we have chosen the point in the phase 
diagram where a fluid mixture with a 1:1 composition 
(that is, xa — 0.50) coexists with a crystal phase that 
is an fee composed of only large particles. We inde- 
pendently calculated the coexistence conditions for this 
point in the phase diagram and we determine the coex- 
isting pressure to be P = 20.1a A /kT , with packing frac- 
tions for the crystal and fluid calculated to be r\ c = 0.61 
and 77/ = 0.51, respectively. 



III. CALCULATION OF INTERFACIAL 
PROFILES 

To monitor changes in the structural or dynamical 
properties across the interface, the system is divided into 
bins along the z-axis, defined as that perpendicular to 
the interfacial plane. Quantities of interest are then cal- 
culated for each bin generating a z-dependent interfacial 
profile for the specific property (density, concentration, 
diffusion, etc.) being measured. The techniques of pro- 
file generation and analysis are similar to those used e. 
lier in the works of Davidchack and Laird on the. singlel 
and binary hard-sphere systems (with a — 0.9)Ej. In this 
section these techniques are summarized, with particular 
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attention to the present calculation. The reader is urged 
to consult the earlier papers if more detail is required. 

In our analysis of the simulations, we employ bins of 
two different resolutions: a coarse scale and a fine scale. 
Coarse scale bins have a width equal to the layer spacing 
of the bulk crystal. This spacing is 0.753(7,4 for [100] and 
0.870cr a for [111]. The fine scale is 1/25 of the coarse 
scale. Fine scale bins reveal in greater detail parame- 
ter variations across the interface, but the coarse scale 
is more useful for observing overall trends in the inter- 
facial profiles. Also for some parameters, such as the 
diffusion constant, only the coarse scale can be used if 
one is to achieve meaningful statistical accuracy. For 
interfacial profiles that exhibit oscillations on the order 
of the lattice spacing, such as density, the conversion 
between the fine-scale profiles and coarse-scale profiles 
to illustrate bulk trends is problematic, since the dis- 
tance between the peaks of such profiles is not necessarily 
constant through the interface. The mismatch between 
the coarse-scala-bins and the peak spacing can lead to 
spurious resultsea if one simply averages over the fine 
scaled bins to create the coarse scaled profile. For such 
profiles we-jemploy a Finite Impulse Response filtering 
procedureCJ to average out the oscillations and reveal 
coarse-grained trends. Details of the specific filtering 
procedure we use can be found in Ref. 24. 

Below is a description on how the various interfacial 
properties were determined. In the definitions, the size 
of the bin is denoted by Az and L x , L y and L z are the 
the dimensions of the simulation box in the x, y and z 
directions, respectively. 

• Pressure: The total pressure profile is defined as 
P(z) = i {P xx (z) + P yy (z) + P zz (z)} , (6) 
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where c indexes the collisions, m is the mass of each 
sphere, < Ef. > is the average kinetic energy per 
sphere, iV c is the number of collisions that occurred 
over the time interval At in the region between 

(c) 

z — Az and z + Az, r k is the kth component of the 
relative distance between the two colliding spheres 

(c) 

and Av k is the kth component of the change in 
velocity for collision c. The first term in Eq. |?] 
represents the ideal gas pressure and the second 
term is the excess part due to sphere interactions. 

• Excess Stress Profiles: The local excess stress is 
calculated from the pressure tensor components. 

S(z)=P zz (z)-±{P xx (z) + P yv (z)} (8) 



In a simulation of an equilibrium interfacial sys- 
tem this quantity should be zero, except in a small 
region at the interface. Improper preparation or 
equilibration of the system often manifests itself in 
the excess of this quantity in the bulk crystal away 
from the interface. As such this quantity is care- 
fully monitored as a measure of the quality of the 
simulation. To smooth out the large oscillations 
in this quantity through the interface, the profile 
is filtered to easily reveal overall trends. (The lo- 
cal excess stress can be integrated with respect to 
z to give the surface excess stress. For a liquid- 
vapor interface the surface excess stress is identical 
to the interfacial free energy, but since the relax- 
ation time for stress in a crystal is generally much 
longer than a typical simulation time, the surface 
excess stress and "f s i canJpe significantly different 
for crystal-melt interfaces!) 

• Density Profiles and Contour Plots: The fine-scale 
density profile for a sphere of type i is determined 
from the number density of that type particle in 
each fine-scale bin. 



Pi{z) 
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where <A^(z)> is the average number of spheres 
of type i in the region between z — Az/2 and z + 
Az/2. To observe overall trends in bulk density 
(or concentration) changes we also produce filtered 
density profiles using our FIR filtering procedure 
discussed above. 

In addition to the z dependent density profiles, it is 
also useful to examine the density variations within 
the x-y planes parallel to the interfacial plane. To 
do this we divide the system into orthorhombic 
subcells with a width in the z direction equal to 
the coarse bin spacing and x and y dimensions of 
0.15(7,4. By counting the average number of par- 
ticles of each type in each subcell and dividing by 
the subcell volume we can produce 2d contour plots 
of the cross-sectional density variation within each 
interfacial plane. 

• Interface location: We determined the location of 
the interface from the orientational order parame- 
ter profile. 



Qn(z) 



jj-^2cos{n9 xy {i,j,k)} 



(10) 



where n is an integer, i,j and k are nearest neigh- 
bor atoms, xy (i,j, k) is the bond angle formed by 
i,j and k projected on the x, y plane, and N z is the 
total number of atoms that form bond angles. The 
average is taken over the number of angles found 
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between z — Az/2 and z + Az/2. The interface in 
the [100] orientation is the point along the z-axis 
where qi(qe for the [111]) is the arithmetic mean 
of the bulk crystal and liquid values. For compar- 
ison, the position of the Gibbs dividing surfaced is 
also calculated. We determine the Gibbs dividing 
surface as the plane along the z-axis such that for 
the 'solute' i, T l = in the equation 

N i /A = p i s z + P i L {L z -z)+T i (11) 

where N l the total number of spheres of type i, A 
is the area of the interface, p l s and p\ are the bulk 
densities, z is the location of the interface assuming 
the length of the simulation box runs from to L z 
and P is the excess particle per unit area of the 
interface. 

• Diffusion coefficient profile: To study the dynam- 
ics across the interface, the diffusion coefficient 
profile is calculated. For a particle of type i, the 
diffusion coefficient is defined as follows 

1 d N ' (z) 

^) = ^) S E(^*)-iN 2 ) (12) 

j=i 

The term in the summation is the mean squared 
displacement over a time interval t — to of a total 
of Ni type i spheres located between z — Az/2 and 
z + Az/2 at time to. 



IV. CONSTRUCTION AND EQUILIBRATION 
OF INTERFACE 

Initially, blocks of crystal and fluid spheres at the cal- 
culated coexistence packing fractions and concentrations 
were prepared separately. As a reference, the z-axis 
is taken to be perpendicular to the interface. The x-y 
planes for both blocks had the same dimensions so that 
they would fit perfectly when put together to construct 
the interface. The plane perpendicular to the interface is 
made as close to square as possible given the geometric 
constraints of the specific interfacial orientation under 
study. This is trivial to achieve with the [100] orienta- 
tion but for [111], the x and y lengths can only be made 
approximately equal. The lengths along z were made 
longer than both those in x and y so that bulk properties 
will be observed between the two interfaces formed. Pe- 
riodic boundary conditions are applied in all directions, 
which results in the two independent crystal-melt inter- 
faces formed along z. The similarity of the two interfaces 
is an important monitor on the quality of the simulation. 
Obviously, if statistically significant differences in struc- 
ture or dynamics exist between the two interfaces, then 
the system has not been properly equilibrated. 



The crystal block with [100] orientation was set up 
with 7776 large spheres. It consisted of 48 crystal layers, 
each layer having 162 spheres. Using the coexistence 
packing fraction T) c — 0.61, the following dimensions 
for the [100] crystal block were used: L x = 13.56(7,4, 
L y = 13.56(7,4 and L z = 36.15a a- Its coexisting fluid 
had 7776 large spheres and 7776 small spheres (15552 
spheres total). The block length is L z — 43.78&A- For 
reasons that will be explained later, this L z gives a 
packing fraction that is slightly higher than that ob- 
tained from the calculated coexistence conditions. For 
the simulation of the [111] interface, the crystal block 
used contained 8190 large spheres, with 45 layers in the 
z direction giving 182 spheres per layer. The crystal 
block dimensions are L x — 13.85(7,4, L y = 12.91(7,4 and 
L z = 39. 13a a- The total number of fluid spheres used 
was also 15552 as in that for the [100] simulation with 
L z = 45.09(7^, again, giving a packing fraction slightly 
higher than that predicted for coexistence. Thus the to- 
tal number of particles in the interface simulations are 
23328 and 23742 for the [100] and [111] interfacial orien- 
tations, respectively. 

Both crystal and fluid blocks are equilibrated sepa- 
rately. The two blocks are then put together but a gap 
equal to a a is left between each of the two crystal-melt 
interfaces formed to ensure that no initial overlap will 
occur at the interfaces. The molecular dynamics simu- 
lation is then started with only the fluid spheres allowed 
to move (the crystal spheres remain fixed). The fluid 
then fills the gaps. To compensate for the decrease in 
the overall bulk density of the fluid phase during this 
step, the fluid blocks are prepared at a packing fraction 
that is slightly higher than the predicted coexistence val- 
ues (as mentioned earlier). In the next step, the crystal 
is equilibrated with the fluid spheres held fixed. At this 
point the interface setup is complete and an equilibration 
run is started with all spheres moving and with initial 
velocities assigned according to a Maxwell distribution. 
In order to efficiently carry out the molecular dynamics 
simulation of .such a large system we use the cell method 
of RappaportEa. 

The stability of a crystal-melt interface in a simulation 
is extremely sensitive to the .assumed coexistence condi- 
tions. In our previous worktja, it was found that the 
predetermined coexistence conditions generally had to 
be modified slightly in order to create a stationary inter- 
face with a zero excess stress in the bulk crystal region. 
This is necessary because a) the coexistence conditions 
are often not known a priori to the accuracy required for 
interface stability and b) the presence of the interface in 
a finite simulation can shift the coexistence equilibrium 
slightly. During our preliminary runs for the current sys- 
tem, using the coexistence conditions as calculated by 
thermodynamic integration of the free energies of sepa- 
rate bulk phases, we found that the resulting interface 
was stable, but yielded a bulk crystal with negative ex- 
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cess stress. Through experimentation, we found that an 
equilibrium interface with zero crystal excess stress was 
possible if the initial fluid packing fraction was increased 
to rjf — 0.52. This had the effect of changing the con- 
centration equilibrium slightly away from a 1:1 mixture 
in the fluid, as discussed below. Now it is in principle 
possible to vary both the initial fluid concentration and 
packing fraction so that the final equilibrium gives pre- 
cisely a 1:1 fluid mixture, however this procedure is quite 
tedious and since our choice of the 1:1 fluid at coexistence 
was arbitrary, the fact that the actual system deviates 
slightly from this concentration is not important for the 
purposes of the current study. 
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FIG. 2: Time evolution of (a) total pressure of the system, (b) 
stress in bulk crystal, (c) fluid densities and (d) location of 
the interfaces, zq. The time unit (mcr^/fceT) 1 / 2 corresponds 
roughly to 18 collisions per particle (cpp). 

To ensure that the system is indeed in equilibrium and 
that the bulk crystal is free of excess stress, we moni- 
tor a variety of properties such as total pressure, bulk 
crystal stress, fluid bulk densities and interfacial loca- 
tion. The results for the [100] interface are shown in 
Fig. |^, which shows that prior to equilibration at about 
t* = t(kT/ma 2 A ) 1/2 = 10000 the crystal grows by about 
3 crystal lattice planes (see Fig. ||d), accompanied by a 
pressure drop from 20.6 to its equilibrium value of 20.1 
a A /ksT (Fig. ^a). In addition, the average excess stress 
in the bulk crystal, initially positive, goes to zero (within 
flucuations) when equilibrium is reached (Fig. [^d). (This 
average excess stress was calculated by averaging S as 
defined above over the middle 28 layers of the bulk crys- 
tal.) 

Initially, the bulk densities of both particle types in the 
fluid are equal, but as the system equilibrates, the bulk 
density of the small particles increases. This increase is 
due to the growth of the crystal (see Fig. ||d). Large 
fluid particles near the crystal freeze, expelling the small 



particles, which are immiscible in the crystal at this pres- 
sure, into the bulk fluid region. Although the bulk fluid 
initially has a large sphere mole fraction of Xa = 0.50, 
the value at equilibrium is somewhat lower (0.46 and 
0.47 for the [100]and [111] interfaces, respectively). The 
equilibrium packing fraction of the bulk fluid slightly re- 
duced from its initial value of rjf — 0.52 to 0.51. Once 
the system is in equilibrium, the interfacial positions are 
stable and the fluctuation in position is less within one 
layer spacing. 

In the preparation of the [100] interface some small 
particles became trapped within some of the interior 
crystal layers as the crystal grew during equilibration. 
Since these were in regions where the diffusion constant 
for the small (and large) particles was found to be zero, 
it cannot be determined whether these particles would 
actually be present in a true equilibrium interface. In 
order to determine the importance of these interstitial 
small particles in stabilizing the interface, we removed 
the particles (about 77 total) from the inner 3 crystal 
layers where they were found. The removal was done at 
t* = 8000 in the equilibration run. Initially the crys- 
tal stress became negative, but quickly returned to zero 
(within fluctuations) as small particles from the bulk dif- 
fused in to reoccupy the removed layer closest to the 
interface (this layer corresponds to layer B in Fig. ||, 
discussed in the next section). The inner two layers 
did not fill in. The interfacial position remained sta- 
ble during this process. The question of true chemical 
equilibrium is always a tricky one in these types of inter- 
face simulationsE^I due to the extremely slow relaxation 
of concentration in the deeper crystal layers. However, 
in this region the concentration of small particles is in 
any event probably quite small and should not affect our 
results significantly (except for perhaps the interfacial 
segregation). As a possible check to this-»roccdure, one 
could use the Widom insertion methodEH to determine 
the excess chemical potential, and thus the solubility, of 
the small particles in the various inner crystal layers, but 
this was not done here. 

The total length of the averaging run after equilibra- 
tion was t* — 4000, which was divided into 40 separate 
blocks of length t* = 100 (corresponds to about 1800 col- 
lisions per particle(cpp)),over which the interfacial pro- 
files were averaged. Since the system contains two inter- 
faces, each block average yields two independent profiles 
(when properly folded about the center of the crystal) 
Thus, each of the profiles reported here represents an 
average of 80 block averages. 

It is important to compare the two independent inter- 
faces produced in a single interface simulation to ensure 
that they are statistically identical. Significant differ- 
ences between the two interfaces are indications of prob- 
lems with the equilibration procedure. As a diagnostic 
we determine the excess stress profile (calculated on the 
fine scale and filtered using the FIR filter described above 
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spheres start developing some order that is similar to the 
large spheres. 
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FIG. 3: Filtered excess stress profiles for the [100] and [111] 
interface orientations. 



and in Ref. 24). These filtered stress profiles are shown 
in Fig. H for both the [100] and [111] orientations - note 
that, the crystal is in the middle of the simulation box. 
The profiles are remarkably symmetric and also show 
that the excess stress is zero within fluctuations in the 
bulk crystal region, ft should be noted that in contrast 
to the case for a liquid-vapor interface, the interfacial 
free energy of a crystal melt interface cannot be deter- 
mined from the integral of the excess stress profile, as the 
relaxation time for stress in the crystal is significantly 
longer than possible simulation timesu, and must be de- 
termined by other means, such as the recently developed 
cleaving wall method!]. The excess stress profiles shown 
in Fig. show a significant negative stress region on the 
crystal side of the interface, indicating that in this region 
the transverse pressure components are greater than the 
pressure component normal to the interfacial plane. The 
precise origin of this unrelaxed crystal stress at the in- 
terface is as yet unknown. 

As mentioned above, the position of the interface is 
determined as the value of z at which the orientational 
order parameter for the large spheres is the arithmetic 
mean of that quantity in the two bulk phases. This 
quantity is a useful measure of interfacial location as 
it is monotonia as a function of z (so that using the 
arithmetic mean makes sense) and can be calculated as 
smooth function without large fluctuations using rela- 
tively short simulation runs. Orientational order param- 
eters 54 and qe, as defined by Eqn. [l^, were determined 
for each particle type. These are shown in Fig. ||. Since 
the crystal phase is made up of pure large spheres and 
we want to see how the ordering of particles is changed 
from bulk crystal to bulk liquid, we determined the in- 
terface location from the parameters calculated for large 
spheres. We also show and gg for the small parti- 
cles and we see that at the interfacial region, the small 
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FIG. 4: Orientational order parameters 54 and qe for the two 
sphere types and for both interfaces studied. 



V. RESULTS FOR THE [100] AND [111] 
INTERFACES 

A. Structure 

The fine scale density profiles for the [100] and [111] 
interfaces are shown in the upper panels of Figures |5| and 
^, respectively. Shown in the lower panels are the cor- 
responding filtered profiles (including the total density 
profile). The distance along the z-axis (in units of the 
large particle diameter, a a) is measured relative to the 
interface center, defined as discussed above by the ori- 
entational order profiles. The vertical dotted lines are 
equally spaced and constructed to correspond with the 
density minima in between the bulk crystal layers. In 
both figures, specific interfacial layers are labeled alpha- 
betically for later reference - layers A and G correspond 
to bulk crystal and liquid, respectively, and layers B-F 
lie within the interfacial region. 

The density profiles for the large particles resemble 
strongly j-those for the single component hard sphere 
interfaced with the periodic oscillations of the bulk crys- 
tal transforming to the uniform density of the fluid over 
about 7-9 lattice layers as the interface is traversed along 
the z-axis. The new feature seen in the present simu- 
lation is the decay of the small particle density over a 
similar distance into the bulk crystal, in which the small 
particle are immiscible. As the small particle density 
decreases into the crystal, it develops oscillations with a 
wavelength closely matching that of the crystal lattice 
spacing. For the [100] interface, the oscillations in the 
small particle density, pb{z) line up in phase with those 
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FIG. 5: Fine-scale (upper panel) and filtered (lower panel) 
density profiles for the [100] orientation. The solid line and 
dashed lines are for the larger (A) and smaller (B) particles, 
respectively. In the lower panel the dotted line shows the 
filtered total density. 




FIG. 6: Fine-scale (upper panel) and filtered (lower panel) 
density profiles for the [111] orientation. The solid line and 
dashed lines are for the larger (A) and smaller (B) particles, 
respectively. In the lower panel the dotted line shows the 
filtered total density. 



so this effect is reminiscent of premelting transitions at 
solid/ vapor interfaces below the bulk melting point, in 
that the presence of a nearby triple point (in this case 
the fcc/NaCl/fluid triple point) manifests itself in the 
presence of the metastable phase (NaCl) at the interface 
between the two coexisting phases (fee and fluid).—. 

As in the single component hard-sphere systemE.il, the 
spacings between the density peaks exhibit some varia- 
tion across the interface - especially for the [100] orienta- 
tion. For each interface, the peak spacing was measured 
by determining the distance between density peaks in the 
fine scale profiles. The resulting peak spacings as func- 
tions of z are shown in Fig. [?]. For the large particles 
the dependence of the spacing on interfacial orientation 
and z is identical to that seen in the single component 
simulations^. The spacing for the [100] lattice increases 
by nearly 20% from the bulk crystal value of 0.76cr a to 
the limiting value of about 0.9(7,4 as the bulk fluid is 
approached. The spacing for the large particles in the 
[111] interface has the same bulk liquid limiting value, 
but since the bulk crystal spacing is very close to this 
limiting value, the variation in spacing across the inter- 
face is quite small. The changes in peak spacing for the 
small particles are quite different for the different ori- 
entations and loosely follow those of the large particle 
- in [100] the small and large particle curves have very 
similar shape, but are shifted by about ua- 



-1 




of the large particle density,/^ (z); whereas, in the [111] 
interface the oscillations are out of phase - the peaks of 
Pb(z) correspond to minima of pa(z). Analysis of the 
atomic positions indicate that this difference is due to 
the fact that in the interfacial region the small particles 
occupy interstitial sites of the large particle fee lattice 
- corresponding to the positions that would be occu- 
pied in an NaCl structure. These preferred positions 
lie in the [100] plane, but lie between the [111] planes 
of the bulk fee lattice. Recall that the NaCl structure 
is the stable structure for this system at high pressure, 



FIG. 7: Peak spacing as determined from maxima of fine- 
scale density profiles for both interfaces studied. 

It is useful to compare these results directly with the 
single component caseE-i In Fig. |^ we plot (upper panel) 
the fine scale density profiles for the [100] orientation of 
both the single component and binary interfaces. The 
single component data was shifted slightly along z to 
make the liquid peaks commensurate. From this plot one 
sees that the presence of the small particles has neglig- 
ble effect on the coexisting liquid density and structure; 
however, the higher pressure for the binary coexistence 
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FIG. 8: Comparison of the binary interface with the ppep 
viously studied hard-sphere single component simulationcj. 
The upper panel shows the [100] fine scale density for both in- 
terfaces. The single component data was shifted along the z- 
axis slightly to maximize the peak overlap in the fluid phase. 
The lower panel shows a comparison with the lattice spacing 
of the [100] interface - for comparison purposes, the data is 
scaled and shifted (vertically) so that all curves go from zero 
in the crystal to unity in the fluid.) 



does give a crystal phase with a higher density (the peaks 
are more closely spaced and more localized). The close 
similarity to the single component system indicates that 
the structure for the large particles is changed very little 
due to the presence of the smaller ones - except for the 
higher density of the crystal. In the lower panel of Fig. ^ 
is shown the peak spacing for the [100] single compo- 
nent and binary interfaces - scaled and shifted so that 
the curves go from zero in the crystal to unity in the 
fluid. The curves for the large particles are qualitatively 
similar, but the change in the single component case is 
less abrupt than that of the binary system. 

A convenient measure of the width of the interfacial 
region is the so-called 10-90 width defined as the dis- 
tance over which an interfacial profile changes from 10% 
to 90% of the higher of the two coexisting bulk values 
relative to the lower bulk value. Such a definition is 
only useful for those interfacial profiles which are mono- 
tonic across the interface, such as a coarse-grained (fil- 
tered) density or diffusion constants. For the filtered 
large particle densities the 10-90 widths are 2. 60^1 for 
the [100] and IAua for the [111] - these are lower by 
about P-8ct,4 than those found for the single component 
systernE-3 which were about 3.3er for the two interfaces. 
From the small particle densities, the widths are larger 
at 3.4<7a and 2>.1oa for the [100] and [111] orientations, 
respectively, The 10-90 region defined by the large par- 
ticles is within that defined by the small particles. The 
larger 10-90 width of the small particle filtered density 



is due to the ability of the small particles to penetrate 
into the first few crystal lattice layers. 
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FIG. 9: Cross-sectional (x-y) density distributions of the large 
spheres for different layers of the [100] interface. 
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FIG. 10: Cross-sectional (x-y) density distributions of the 
small spheres for different layers of the [100] interface. 

To get a more detailed picture of the transition from 
crystal-like to fluid-like structure as the interface is tra- 
versed it is useful to examine the density distributions 
within x-y cross-sectional planes parallel to the interface. 
(The reported distributions are averages taken over 1800 
epp - details of their calculation can be found in the 
previous section.) Figures [)] and [lO] respectively show 
the x-y large and small particle density distributions for 
the [100] interface orientation as greyscale contour plots. 
The layer labels A-G correspond to those shown in Fig. [|. 
Fig. H shows that this transition from crystal to fluid oc- 
curs over about three layers (C,D and E) for the [100] 
interface and that these transition layers are not uniform, 
but consist of coexisting solid- and liquid-like regions, as 
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was seen in the single-component simulationscJ. Layer 
B, although fully crystalline, does possess two vacancy 
defects at points (-3.3, -3.0) and (-1.0,5.3). The [100] 
contour plots for the small particle density are quite in- 
teresting. There is considerable density in layer B where 
the small particles are present in two types of positions 
- in the 'NaCl' interstitial positions and in the positions 
corresponding to the vacancies of the large particle crys- 
tal lattice found in layer B. The interstitial positions are 
occupied by single small particles, but each vacancy is 
filled with severaLsmall particles. In the single compo- 
nent simulationsE-3 vacancy nucleation at the interface 
was also seen, in that case the vacancies once formed 
were highly mobile, migrating into the bulk via a hop- 
ping mechanism. In the present simulations, however, 
once the vacancies are formed in the large particle lat- 
tice, they are quickly filled with some number of small 
particles, which appears to immobilize the defect by sup- 
pressing the hopping mechanism - however the evidence 
for this is anecdotal, as the number of such vacancies is 
too small to gather meaningful statistics. 

To estimate the degree of interfacial segregation, the 
Gibbs dividing surface for both interfacial orientations 
was determined according to Eq. [ll] and found it to be 
close to the interface location determined from the orien- 
tational order parameters. The surface is at z = — Q.buA 
for [100] and at z = ~0.9a A for [111]. At these divid- 
ing surfaces, the excess density of solute (here defined 
as component B) was found to be negligible - indicating 
minimal interfacial segregation. Of course, for such in- 
terfacial simulations, the question of complete chemical 
equilibrium is generally problematic, as discussed in the 
previous section; however, we are confident that the con- 
centrations of each particle type from interfacial layer B 
out to the bulk fluid are in chemical equilibrium (since 
diffusion is non-negligible there) and that the equilib- 
rium concentrations of small particles in layers deeper 
into the crystal are probably quite small and will not 
significantly affect the results presented here. 
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FIG. 11: (a) Diffusion coefficient profile for the [100] inter- 
face, (b) Scaled diffusion coefficients. 
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FIG. 12: (a) Diffusion coefficient profile for the [111] inter- 
face, (b) Scaled diffusion coefficients. 



B. Dynamics 

We study the dynamics across the interface by measur- 
ing diffusion coefficients in the coarse-scaled bins. The 
diffusion profiles for the [1001 and [111] interfaces are 
shown in Figures pd](a) and |l2|(a), respectively. 

The limiting bulk diffusion coefficient is 
0.012(fcsTcr^/m) 1 / 2 for the large spheres and 
0.050(/c B T(T^/m) 1 / 2 for the small particles, inde- 
pendent of the crystal orientation, as expected. When 
the three Cartesian components of the total diffusion 
coefficient are separately determined, it is found that 
diffusion is isotropic throughout the interfacial region. 

The larger value of the small particle diffusion con- 
stant makes it difficult to compare the diffusion constants 



of the two components so we also plot for each, the ra- 
tio diffusion constant to the average fluid bulk value in 
Figures |ll](b) and |l|(b). Here we find the interesting 
result that the two curves (for both crystal orientations) 
are similar in shape, but shifted relative to one another 
by more than la a- As the interface is traversed from 
fluid to crystal, the diffusion constant for the large par- 
ticle goes effectively to zero near z = 0, but the small 
particles still have significant mobility. In this region, the 
large particles have become "locked in" to their crystal 
lattice sites, but the small particles can still move about 
- primarily by hopping between interstitial sites. 

The 10-90 widths from the diffusion coefficient profiles 
for both orientation and particle types are about 3<t a- 
But because the diffusion profiles are shifted, the 10-90 
widths do not define the same region. If contributions 
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from both particle types are considered, the widths are 
4.5a A for the [100] and 3.9a A for the [111] interface. The 
center of these interfacial regions are shifted by about 
la a to the fluid side compared to the interfacial regions 
defined by the density profiles. To illustrate this more 
clearly, we show in Fig. ^ all of the order parameter 
profiles (orientation, diffusion and density) for the [100] 
interface, scaled in such a way that they go from unity 
in the crystal phase to zero in the liquid (for example for 
the diffusion constants we plot 1 — D(z)/Df). 
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FIG. 13: Diffusion, orientation and filtered density order pa- 
rameter profiles for the [100] interface - all scaled such that 
they go from unity in the crystal to zero in the fluid phase. 

The 10-90 regions for the diffusion constants are offset 
(toward the liquid side) from those for the filtered density 
profiles so the interfacial region is wider than any single 
structural or dynamical quantity would indicate. If one 
considers the interfacial region as the union of the 10-90 
regions for the separate profiles, then the width of the 
interfacial region is 4.8(7,4, greater than that calculated 
from densities or diffusion coefficients alone. 



VI. SUMMARY 

We have performed a series of molecular-dynamics 
simulations to study the crystal-melt interface of a bi- 
nary hard-sphere system with diameter ratio 0.414. Pre- 
vious simulation studies on two-component crystal-melt 
interfaces have focused on equilibrium between a. fluid 
mixture and a substitutional^ disordered crystaHE3, 
but here we have examined the interface between a fluid 
mixture (approximately equimolar in concentration) and 
a coexisting smgZe-component fee crystal comprised of 
large particles - in which the small particles are im- 
miscible. Such a coexistence occurs at relatively low 
pressures in the phase diagram for this diameter ratio 
- at higher pressures the fluid coexists with a 1:1 or- 



dered crystal with an 'NaCl' structure. At a pressure of 
P = 20.la A /kT the two phases coexist at the following 
packing fractions: r/ c — 0.61 and rjf = 0.51. 

Some of the principal results of this study are as fol- 
lows: 

• The interfacial density profiles of the large parti- 
cles is very similar to that of the single-reomponent 
hard-sphere system previously studiedEa, indicat- 
ing that the presence of the small particle has no 
significant effect on the interfacial structure of the 
large particle, except for a compression of the crys- 
tal lattice due to the higher pressure. In particular 
the variation of the spacing between the large par- 
ticle density peaks is very similar to that found in 
the single component studies. 

• Within the regions of the interface in which the 
large particles are largely confined to fee lattice 
sites, the small particles occupy either vacancy 
sites in the fee lattice or 'NaCl' interstitial sites. 
The interstitial sites are single occupied whereas 
the vacancy sites are found to be occupied by sev- 
eral small particles. The presence of the small par- 
ticles greatly suppresses the mobility of the fee va- 
cancies relative to those previously Bated in single- 
component hard-sphere simulations^. 

• There does not appear to be significant solute 
(small particle) segregation at the interface. 

• The diffusion profiles of the small and large par- 
ticles are similar in width (about 3 <ja), but are 
shifted relative to one another by about 1 a a along 
the interface normal (z axis). As one traverses the 
interface from bulk fluid to bulk crystal, the dif- 
fusion constant goes to zero for the large particles 
in a region in which there is still significant small 
particle mobility. The picture in this region is of 
large particles localized at fee lattice sites, with 
the small particles still diffusing between intersti- 
tial site within the lattice of large particles. 

• As was .found in previous hard-sphere interface 
studicsE3t3 the total width of the interfacial re- 
gion is greater than the width determined by any 
single interfacial profile (such as diffusion or den- 
sity) as the profiles for the individual quantities can 
be significantly shifted from one another. Specifi- 
cally we see that as one moves from the crystal into 
the fluid the bulk density relaxes first to liquid-like 
values before significant mobility (diffusion) is ob- 
served. Considering both structural and dynamic 
properties, the interfacial (10-90) width is 4.8cr^. 
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